# Hello, world!
#
# This is an example function named 'hello'
# which prints 'Hello, world!'.
#
# You can learn more about package authoring with RStudio at:
#
#   http://r-pkgs.had.co.nz/
#
# Some useful keyboard shortcuts for package authoring:
#
#   Install Package:           'Ctrl + Shift + B'
#   Check Package:             'Ctrl + Shift + E'
#   Test Package:              'Ctrl + Shift + T'
#' Title kegg_enrich_local_compareCluster
#'
#' @param gene_list              a gene_list for enrichment, the gene list of length range from 1-max
#'
#' @return
#' @export
#' @importFrom data.table fread
#' @examples   x<-kegg_enrich_local_compareCluster(species,gene_list[[1]],outdir)
hello <- function() {
  print("Hello, world!")
}

#' Title mkdirs mkdir new dirs without overwrite
#' @param outdir  which dir for create new dirs
#' @param fp      The name for new dirs
#'
#' @return
#' @export
#' @examples  mkdirs(outdir,new_dirs_name)
mkdirs <- function(outdir,fp) {
  if(!file.exists(file.path(outdir,fp))) {
    #		mkdirs(dirname(fp))
    dir.create(file.path(outdir,fp))
  }else{
    print(paste(fp,"Dir already exists!",sep="     "))
    #unlink(file.path(outdir,fp), recursive=TRUE)
    #dir.create(file.path(outdir,fp))
  }
}
#' Title mkdirs mkdir new dirs with overwrite
#' @param outdir  which dir for create new dirs
#' @param fp      The name for new dirs
#'
#' @return
#' @export
#' @examples  mkdirs_over(outdir,new_dirs_name)
mkdirs_over <- function(outdir,fp) {
  if(!file.exists(file.path(outdir,fp))) {
    #		mkdirs(dirname(fp))
    dir.create(file.path(outdir,fp))
  }else{
    print(paste(fp,"Dir already exists!",sep="     "))
    unlink(file.path(outdir,fp), recursive=TRUE)
    dir.create(file.path(outdir,fp))
  }
}
#' Title cell_vdj cell with vdj filter
#' @param cell  which dir for create new dirs
#' @param chain chain of XCR default[chain=c('TRB','TRA')],
#'
#' @return data.frame
#' @export
#' @importFrom dplyr arrange
#' @examples  cell_vdj<(cell_vdj_data.frame,chain=c('TRB','TRA'))
#'
cell_vdj<-function(cell,chain=c('TRB','TRA')){
  cell_vdj<-data.frame(raw_clonotype_id=0,barcode=0,trb_v_gene=0,trb_d_gene=0,trb_j_gene=0,trb_c_gene=0,trb_cdr3=0,tra_v_gene=0,tra_d_gene=0,tra_j_gene=0,tra_c_gene=0,tra_cdr3=0)
  for(cell_barcode in unique(cell$barcode)){
    trb<-c()
    tra<-c()
    #cell_barcode<-"ACACTGAAGTGACATA-1"
    tmp2<-cell[cell$barcode==cell_barcode,]
    tmp2<-dplyr::arrange(tmp2,desc(umis),d_gene)
    if( dim(tmp2[tmp2$chain==chain[1],])[1]>=1 ){
      tmp<-tmp2[tmp2$chain==chain[1],]
      trb<-c(as.character(tmp[1,c('v_gene')]),as.character(tmp[1,c('d_gene')]),as.character(tmp[1,c('j_gene')]),as.character(tmp[1,c('c_gene')]),as.character(paste(chain[1],as.character(tmp[1,c('cdr3')]),sep=":")))
    }else if(dim(tmp2[tmp2$chain==chain[1] ,])[1]==0 ){
      trb<-c('None','None','None','None','None')
      # }		else if(dim(tmp2[tmp2$chain==chain[1]& tmp2$d_gene=="None",])[1]>=1){
      # tmp<-tmp2[tmp2$chain==chain[1],]
      # trb<-c(as.character(tmp[1,c('v_gene')]),as.character(tmp[1,c('d_gene')]),as.character(tmp[1,c('j_gene')]),as.character(tmp[1,c('c_gene')]),as.character(paste(chain[1],as.character(tmp[1,c('cdr3')]),sep=":")))
    }else {
      print(tmp2)
    }
    #
    if(dim(tmp2[tmp2$chain==chain[2],])[1]>=1){
      tmp<-tmp2[tmp2$chain==chain[2],]
      tra<-c(as.character(tmp[1,c('v_gene')]),as.character(tmp[1,c('d_gene')]),as.character(tmp[1,c('j_gene')]),as.character(tmp[1,c('c_gene')]),as.character(paste(chain[2],as.character(tmp[1,c('cdr3')]),sep=":")))
    }else if(dim(tmp2[tmp2$chain==chain[2] ,])[1]==0){
      tra<-c('None','None','None','None','None')
      # }else if( dim(tmp2[tmp2$chain==chain[2],])[1]==1 ){
      # tmp<-tmp2[tmp2$chain==chain[2],]
      # tra<-c(as.character(tmp[1,c('v_gene')]),as.character(tmp[1,c('d_gene')]),as.character(tmp[1,c('j_gene')]),as.character(tmp[1,c('c_gene')]),as.character(paste(chain[2],as.character(tmp[1,c('cdr3')]),sep=":")))
    }else {
      print(tmp2)
    }
    cell_vdj<-rbind(cell_vdj,a=c(as.character(tmp$raw_clonotype_id[1]),as.character(cell_barcode),trb,tra))
  }
  cell_vdj<-cell_vdj[2:nrow(cell_vdj),]
  return(cell_vdj)
}
#cell_vdj2<-cell_vdj(cell,chain=c('TRB','TRA'))
#' Title trc merge filtered_contig_annotations and clone file with filter
#' @param tcr_folder   which dir for create new dirs
#' @param clone_folder chain of XCR default[chain=c(\'TRB\',\'TRA\')],
#' @param sample       sample name for tcr or bcr
#' @param trim         trim the cellbarcode suffix,default \\-1
#' @param chain        chain of xcr
#' @return data.frame
#' @export
#' @importFrom dplyr arrange
#' @examples  trc (tcr_folder,i,chain=chain)
#'
trc <- function(tcr_folder,sample,trim="\\-1",chain=c('TRB','TRA')){
  #tcr_folder<-"E:/姚盟成/正在进行的任务/免疫VDJ-immunarch/outs/"
  tcr <- read.csv(paste(tcr_folder,paste(sample,".filtered_contig_annotations.csv",sep=""), sep="/"))
  tcr$barcode <- gsub(trim, "", tcr$barcode)
  tcr$barcode <- paste(tcr$barcode,sample,sep=".")
  dim(tcr)
  #[1] 5091   18
  #Also, it may be worthwhile to remove/ignore columns from your VDJ data that are FALSE or None in the "productive" column, as you'll end up with a lot of extra noise otherwise. Actually, it might be best to just get the clonotype id for each cell, then get the gene information from the consensus_annotations.csv file.
  tcr<-tcr[tcr$productive=="True" & tcr$is_cell=="True"& tcr$chain!="Multi",]
  dim(tcr)
  #cell<-cell[cell$is_cell=="True" & cell$high_confidence=="True" & cell$productive=="True",]
  #
  #tcr <- tcr[!duplicated(tcr$barcode), ] remove this comand
  tcr <-cell_vdj(tcr,chain=chain)
  dim(tcr)
  print(head(tcr))
  # Only keep the barcode and clonotype columns.
  # We'll get additional clonotype info from the clonotype table.
  #tcr <- tcr[,c("barcode", "raw_clonotype_id")]
  #tcr <- tcr[,c("barcode", "raw_clonotype_id","chain","v_gene","d_gene","j_gene")]
  names(tcr)[names(tcr) == "raw_clonotype_id"] <- "clonotype_id"
  # Clonotype-centric info.
  clono <- read.csv(paste(tcr_folder,paste(sample,".clonotypes.csv",sep=""), sep="/"))
  # Slap the AA sequences onto our original table by clonotype_id.
  tcr <- merge(tcr, clono[, c("clonotype_id", "cdr3s_aa",'cdr3s_nt')],by="clonotype_id")
  # Reorder so barcodes are first column and set them as rownames.
  #tcr <- tcr[, c(2,1,3:7)]
  #tcr <- tcr[, c(2,1,3)]
  rownames(tcr) <- tcr$barcode
  tcr$barcode<-NULL
  return (tcr)
}
#' Title all_samples_tcr  merge all samples tcr data.frame
#' @param all_samples   all samples for tcr merge
#' @param tcr_folder    tcr file dirs
#' @param chain         chain of xcr
#' @return              data.frame
#' @export
#' @importFrom          dplyr arrange
#' @examples  all_samples_tcr (all_samples,tcr_folder,chain=c('TRB','TRA'))
#'
all_samples_tcr<-function(all_samples,tcr_folder,chain=c('TRB','TRA')){
  all_tcr<-data.frame()
  for (i in all_samples){
    print(i)
    yao<-trc (tcr_folder,i,chain=chain)
    yao$Sample<-rep(i,dim(yao)[1])
    all_tcr<-rbind(all_tcr,yao)
  }
  return(all_tcr)
}
#' Title integration    single cell RNA integration with different samples
#' @param object_list   all samples seurat list
#' @param cca_dims      reduction dims of cca in integration
#' @param pca_dims      reduction dims of pca in integration
#' @param RunPCA_npcs   RunPCA_npcs
#' @return              s4 for seurat
#' @export
#' @importFrom   Seurat FindIntegrationAnchors
#' @importFrom   Seurat IntegrateData
#' @importFrom   Seurat ScaleData
#' @importFrom   Seurat RunPCA
#' @importFrom   Seurat DefaultAssay
#' @examples  integration(object_list)
#'
integration <- function(object_list,cca_dims=20,pca_dims=20,RunPCA_npcs=30){
  all_samples<-names(object_list)
  if(length(all_samples)>1){
    immune.anchors <- FindIntegrationAnchors(object.list = object_list, dims = 1:as.numeric(cca_dims))
    immune.combined <- IntegrateData(anchorset = immune.anchors, dims = 1:as.numeric(pca_dims))
    DefaultAssay(immune.combined) <- "integrated"
  }else{
    print(paste('You just have only one sample: hahahahah',all_samples,sep=' '))
    immune.combined=object_list[[1]]
  }
  immune.combined <- ScaleData(immune.combined, verbose = FALSE)
  immune.combined <- RunPCA(immune.combined, npcs = as.numeric(RunPCA_npcs), verbose = FALSE)

  return (immune.combined)
}

#' Title sample2SeuratObject_list import single cell RNA expression file as seurat
#' @param indir                   which dirs import with expression file
#' @param oudir                   output dirs for files
#' @param samplenames             all samples for import
#' @return                        s4 for seurat list
#' @export
#' @examples                    sample2SeuratObject_list(indir,oudir,samplenames)
#'
sample2SeuratObject_list <- function(indir,oudir,samplenames){
  print("start read")
  SeuratObject_list <- mapply(read_onesample2SeuratObject, indir, oudir,samplenames)
  print (length(SeuratObject_list))
  print("end read")
  names(SeuratObject_list) <- samplenames
  return(SeuratObject_list)
}
#' Title read_onesample2SeuratObject         Import single sample expresssion file as seurat
#' @param indir                              Which dirs import with expression file
#' @param oudir                              Output dirs for files
#' @param samplenames                        All samples for import
#' @param min.cells                          Min cells for gene detected
#' @param mitoName                           Prefix of mit name
#' @param HB                                 HB threshold
#' @param mt.percent                         Threshold of mt.percent
#' @param min_nFeature_RNA                   Min genes for cell detected
#' @param nfeatures_FindVariableFeatures     VariableFeatures for ananlysis
#' @param FindVariableFeatures_method        Find VariableFeatures method
#'
#' @return                                   s4 for seurat list
#' @export
#' @importFrom   Seurat Read10X
#' @importFrom   Seurat CreateSeuratObject
#' @importFrom   Seurat PercentageFeatureSet
#' @importFrom   Seurat VlnPlot
#' @importFrom   Seurat NormalizeData
#' @importFrom   Seurat FindVariableFeatures
#' @importFrom   Seurat FeatureScatter
#' @importFrom   Seurat VariableFeatures
#' @importFrom   Seurat VariableFeaturePlot
#' @importFrom   Seurat CombinePlots
#' @importFrom   Seurat LabelPoints
#' @importFrom   ggplot2  theme
#' @examples     read_onesample2SeuratObject(indir,outdir,sample_names)
#'
read_onesample2SeuratObject <- function(indir,outdir,sample_names,min.cells=5,mitoName="MT",HB=100,mt.percent=10,min_nFeature_RNA=500,max_nFeature_RNA=2500,nfeatures_FindVariableFeatures=2000,FindVariableFeatures_method="vst"){
  print(paste('Start read rawdata, the sample is :',sample_names,sep=' '))
  sample_indir<- file.path(indir, sample_names,'/outs/filtered_gene_bc_matrices/ref/')
  #expression_matrix <- Read10X(data.dir = sample_indir)
  sample_name.data<-data.frame()
  if (file.exists(sample_indir)){
    sample_name.data <- Read10X(data.dir = sample_indir)}
  else{
    sample_name.data <- read.table(file.path(indir,paste(sample_names,'.csv',sep='')), sep = ",",header=TRUE,row.names = 1)
  }
  colnames(sample_name.data)<-paste(colnames(sample_name.data),sample_names,sep=".")
  SeuratObject <- CreateSeuratObject(counts = sample_name.data, project = sample_names, min.cells = min.cells, min.features = 200)
  SeuratObject$stim <- sample_names
  SeuratObject[["percent.mt"]] <- PercentageFeatureSet(SeuratObject, pattern = paste("^", mitoName,"-",sep=""))
  if (sum(SeuratObject@meta.data$percent.mt)==0){SeuratObject[["percent.mt"]] <- PercentageFeatureSet(SeuratObject, pattern = paste("^", 'mt',"-",sep=""))}
  p<-theme(panel.grid=element_blank(),axis.text.x = element_text(color="black",face="plain",size=20), axis.text.y = element_text(color="black",face="plain",size=20), axis.title.x = element_text(face="plain", color="black",size=25),axis.title.y = element_text(face="plain", color="black",size=25),legend.text=element_text(color="black",face="plain",size=15),legend.title=element_text(color="black",face="plain",size=15),strip.text=element_text(face="plain",size=0),title=element_text(size=20),legend.position="top")
  pdf(paste(sample_names,'qc.pdf',sep='_'),w=12,h=8)
  p0<-VlnPlot(SeuratObject, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size =-1)+ggtitle("Before QC(percent.mt,min_nFeature_RNA)")+theme(plot.title = element_text(hjust = 0.5))+p
  A<<-min_nFeature_RNA
  B<<-max_nFeature_RNA
  mt<<-mt.percent
  SeuratObject <- subset(SeuratObject, subset = nFeature_RNA > A & nFeature_RNA < B & percent.mt < mt)
  SeuratObject <-NormalizeData(object = SeuratObject, verbose = FALSE)
  SeuratObject<- FindVariableFeatures(object = SeuratObject, selection.method = FindVariableFeatures_method, nfeatures =nfeatures_FindVariableFeatures)
  setwd(outdir)
  #####质控结果pbmc
  p1<-VlnPlot(SeuratObject, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size =-1)+ggtitle("After QC(percent.mt,min_nFeature_RNA)")+theme(plot.title = element_text(hjust = 0.5))+p
  plot1 <- FeatureScatter(SeuratObject, feature1 = "nCount_RNA", feature2 = "percent.mt")
  plot2 <- FeatureScatter(SeuratObject, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
  p2<-CombinePlots(plots = list(plot1, plot2))
  top10 <- head(VariableFeatures(SeuratObject), 10)
  plot1 <- VariableFeaturePlot(SeuratObject)
  plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE,xnudge = 0, ynudge = 0)
  p3<-CombinePlots(plots = list(plot1, plot2))
  print(p0)
  print(p1)
  print(p2)
  print(p3)
  dev.off()
  #####
  print(paste('Finished read rawdata, the sample is :',sample_names,sep=' '))
  return(SeuratObject)
}
#' Title read_onesample2SeuratObject         Import single sample expresssion file as seurat
#' @param immune.combined                    Seurat list
#' @param outdir                             Output dirs for files
#' @param pref                               Output file name prefix
#' @param w_h                                w and h for pdf file
#'
#' @return                                   s4 for seurat list
#' @export
#' @importFrom   Seurat JackStraw
#' @importFrom   Seurat ScoreJackStraw
#' @importFrom   Seurat DimPlot
#' @importFrom   Seurat VlnPlot
#' @importFrom   Seurat DimHeatmap
#' @importFrom   Seurat VizDimLoadings
#' @importFrom   Seurat JackStrawPlot
#' @importFrom   Seurat ElbowPlot
#' @importFrom   ggplot2  theme
#' @examples     qc_pca_plot(immune.combined)
#'
qc_pca_plot<-function(immune.combined,outdir=getwd(),pref='10x',w_h=c(12,8)){
  # plot pca
  pbmc <- JackStraw(immune.combined, num.replicate = 100)
  pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
  pdf(paste(outdir,paste(pref,"pca.qc.pdf",sep='_'),sep='/'),w=w_h[1],h=w_h[2])
  p1<-DimPlot(immune.combined, reduction = "pca")
  # plot pca
  #    pdf(paste(outdir,paste(pref,"pca_heatmap.pdf",sep='_'),sep='/'),w=w_h[1],h=w_h[2])
  p2<-DimHeatmap(immune.combined, dims = 1:15, cells = 500, balanced = TRUE)
  #    dev.off()
  # plot pca
  #    pdf(paste(outdir,paste(pref,"pca_ElbowPlot.pdf",sep='_'),sep='/'),w=w_h[1],h=w_h[2])
  p3<-VizDimLoadings(immune.combined, dims = 1:2, reduction = "pca")
  p4<-JackStrawPlot(pbmc, dims = 1:15)
  p5<-ElbowPlot(immune.combined)
  print(p1)
  print(p2)
  print(p3)
  print(p4)
  print(p5)
  dev.off()
}
#' Title                                     cluster for single cell RNA with viarable genes
#' @param object_list                        Seurat list
#' @param outdir                             Output dirs for files
#' @param prefix                               Output file name prefix
#' @param w_h                                w and h for pdf file
#'
#' @return                                   s4 for seurat list
#' @export
#' @importFrom   Seurat  RunUMAP
#' @importFrom   Seurat  FindNeighbors
#' @importFrom   Seurat  FindClusters
#' @importFrom   Seurat  Idents
#' @importFrom   plyr    mapvalues
#' @importFrom   Seurat  DimPlot
#' @importFrom   cowplot plot_grid
#' @importFrom   Seurat  ElbowPlot
#' @importFrom   Seurat  DefaultAssay
#' @importFrom   ggplot2  theme
#' @examples     clusters(object_list,outdir)
#'
clusters<-function(object_list,outdir,prefix="clusters",w_h=c(24,8),selection.method = "vst",nfeatures = 2000,anchors_ims=20,group="stim",resolution=0.5,cca_dims=20,pca_dims=20,RunPCA_npcs=30){
  ###进行降维聚类 批次校正等分析
  immune.combined<-integration(object_list,cca_dims=cca_dims,pca_dims=pca_dims,RunPCA_npcs=RunPCA_npcs)
  # UMAP and Clustering
  immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:20)
  immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:20)
  immune.combined <- FindClusters(immune.combined, resolution = as.numeric(resolution))
  #改聚类号，从1开始
  current.cluster.ids <- levels(Idents(immune.combined))
  new.cluster.ids <- as.numeric(current.cluster.ids)+1
  Idents(immune.combined) <- plyr::mapvalues(x = Idents(immune.combined), from = current.cluster.ids, to = new.cluster.ids)
  immune.combined@meta.data$seurat_clusters<-Idents(immune.combined)[rownames(immune.combined@meta.data)]
  print(head(immune.combined@meta.data))
  #画图
  library(cowplot)
  p_tmp<<-theme(panel.grid=element_blank(), legend.background = element_rect(colour = NA),
                legend.title = element_blank(),legend.text =  element_text(color="black",size=30),
                axis.text.x = element_text(color="black",size=30),
                axis.text.y = element_text(color="black",size=30),
                axis.title.x = element_text(face="plain", color="black",size=30),
                axis.title.y = element_text(face="plain", color="black",size=30))
  pdf(paste(outdir,paste(prefix,"umap_cluster_samples.pdf",sep="_"),sep="/"),w=w_h[1],h=w_h[2])
  p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim",pt.size = 1)+p_tmp
  p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE,label.size = 8,pt.size = 1)+p_tmp
  p3<-plot_grid(p1, p2)
  print(p3)
  dev.off()
  return (immune.combined)
}

#' Title  combine_gene_xcr                  combine gene and xcr
#' @param object_list                       Seurat object
#' @param outdir                            tcr_data frame
#'
#' @return                                   s4 for seurat object
#' @export
#' @importFrom   Seurat  AddMetaData
#' @importFrom   Seurat  FindNeighbors
#' @importFrom   Seurat  FindClusters
#' @importFrom   Seurat  Idents
#' @importFrom   plyr    mapvalues
#' @importFrom   Seurat  DimPlot
#' @importFrom   cowplot plot_grid
#' @importFrom   Seurat  ElbowPlot
#' @examples     combine_gene_xcr(seurat_object,tcr_data)
#'
combine_gene_xcr<-function(combined,tcr_data){
  combined<-AddMetaData(combined, metadata = tcr_data)
  combined@meta.data$TRA<-paste(combined@meta.data$tra_v_gene,combined@meta.data$tra_j_gene,sep="_")
  combined@meta.data$TRB<-paste(combined@meta.data$trb_v_gene,combined@meta.data$trb_j_gene,sep="_")
  combined@meta.data$TCR<-as.character(combined@meta.data$cdr3s_aa)
  combined@meta.data$TCR[is.na(combined@meta.data$TCR=="")]<-"No TCR"
  combined@meta.data$TCR[combined@meta.data$TCR!="No TCR"]<-"TCR"
  combined@meta.data$cdr3_aa<-paste(combined@meta.data$tra_cdr3,combined@meta.data$trb_cdr3,sep=";")
  combined@meta.data$cdr3<-combined@meta.data$cdr3_aa
  combined@meta.data[grep("NA|None",combined@meta.data$cdr3_aa),]$cdr3<-"None"
  combined@meta.data[-grep("NA|None",combined@meta.data$cdr3_aa),]$cdr3<-combined@meta.data[-grep("NA|None",combined@meta.data$cdr3_aa),]$cdr3_aa
  combined@meta.data$vdj<-paste(combined@meta.data$TRA,combined@meta.data$TRB,sep=':')
  return(combined)
}

#' Title                                    umap plot for VDJ and seurat_clusters
#' @param object_list                       Seurat object
#' @param outfilename                       Output file name
#'
#' @return                                   s4 for seurat object
#' @export
#' @importFrom   Seurat  AddMetaData
#' @importFrom   Seurat  DimPlot
#' @importFrom   Seurat  FindClusters
#' @importFrom   Seurat  Idents
#' @importFrom   plyr    mapvalues
#' @importFrom   Seurat  DimPlot
#' @importFrom   cowplot plot_grid
#' @importFrom   Seurat  ElbowPlot
#' @importFrom   ggplot2  theme
#' @examples     umap_plot(seurat_object,outfilename)
#'
umap_plot<-function(immune.combined_umap,outfilename='gene_vdj_umap_cluster_TCR.pdf',group.by = "TCR",cols=c("blue","red"),w_h=c(24,10)){
  p_tmp<-theme(panel.grid=element_blank(), legend.background = element_rect(colour = NA),
               legend.title = element_blank(),legend.text =  element_text(color="black",size=30),
               axis.text.x = element_text(color="black",size=30),
               axis.text.y = element_text(color="black",size=30),
               axis.title.x = element_text(face="plain", color="black",size=30),
               axis.title.y = element_text(face="plain", color="black",size=30))
  pdf(outfilename,w=w_h[1],h=w_h[2])
  # Visualization
  p1 <- DimPlot(immune.combined_umap, reduction = "umap", group.by = group.by,pt.size = 1,cols=cols)+p_tmp
  p2 <- DimPlot(immune.combined_umap, reduction = "umap", label = TRUE,label.size = 8,pt.size = 1)+p_tmp
  p3<-plot_grid(p1, p2)
  print(p3)
  dev.off()
}

#' Title  top10_cdr3_plot                   top10_cdr3 plot circ
#' @param top10_cdr3                        cdr3 sequnce for plot
#' @param rds                               seurat object
#' @param outdir                            Outdir for outfile, it will create new dirs with name of top_cdr3_plot
#' @param prefix                            Output file name prefix
#'
#' @return
#' @export
#' @importFrom   Seurat  AddMetaData
#' @importFrom   Seurat  DimPlot
#' @importFrom   Seurat  FindClusters
#' @importFrom   Seurat  Idents
#' @importFrom   plyr   mapvalues
#' @importFrom   Seurat  DimPlot
#' @importFrom   cowplot plot_grid
#' @importFrom   Seurat  ElbowPlot
#' @importFrom   ggplot2  theme
#' @examples     top10_cdr3_plot(top10_cdr3,rds,outdir,prefix)
#'
top10_cdr3_plot <- function(top10_cdr3,rds,outdir,prefix){
  mkdirs(outdir,'top_cdr3_plot')
  setwd(paste(outdir,'top_cdr3_plot',sep='/'))
  outdir<-paste(outdir,'top_cdr3_plot',sep='/')
  for (i in top10_cdr3){
    tmp_rds<-rds #small_immune.combined_umap
    tmp_rds@meta.data$TCR<-rep("NOT IN",nrow(tmp_rds@meta.data))
    tmp_rds@meta.data$TCR<-as.character(tmp_rds@meta.data$TCR)
    tmp_rds@meta.data$cdr3s_aa<-as.vector(tmp_rds@meta.data$cdr3s_aa) #which(tmp_rds@meta.data$cdr3s_aa==i)
    tmp_rds@meta.data$TCR[which(tmp_rds@meta.data$cdr3s_aa==i)]<-rep(paste("top",which(top10_cdr3==i),sep=""),length(which(tmp_rds@meta.data$cdr3s_aa==i)))
    print(unique(tmp_rds@meta.data$TCR))

    p_tmp<-theme(panel.grid=element_blank(), legend.background = element_rect(colour = NA),
                 legend.title = element_blank(),legend.text =  element_text(color="black",size=30),legend.position="right",
                 axis.text.x = element_text(color="black",size=30),
                 axis.text.y = element_text(color="black",size=30),
                 axis.title.x = element_text(face="plain", color="black",size=30),
                 axis.title.y = element_text(face="plain", color="black",size=30))
    pdf(paste(outdir,paste(prefix,which(top10_cdr3==i),"umap_cluster_TCR.pdf",sep='_'),sep='/'),w=24,h=10)
    # Visualization
    p1 <- DimPlot(tmp_rds, reduction = "umap", group.by = "TCR",pt.size = 1,cols=c("grey","red"))+p_tmp
    p2 <- DimPlot(tmp_rds, reduction = "umap", label = TRUE,label.size = 8,pt.size = 1)+p_tmp
    p3<-plot_grid(p1, p2)
    print(p3)
    dev.off()
    print(paste(outdir,paste(prefix,which(top10_cdr3==i),"umap_cluster_TCR.pdf",sep='_'),sep='/'))
    if (dim(tmp_rds@meta.data[tmp_rds@meta.data$TCR!="NOT IN",c('seurat_clusters','TCR')])[1]=="0"){
      print(paste("There is no cdr3 in this sample, the cdr3 is :",i,sep=""))
      next()
    }else{
      print(paste("There is cdr3 in this sample.................................. the cdr3 is :",i,sep=""))
    }
    tmp_a<-as.data.frame(table(tmp_rds@meta.data[tmp_rds@meta.data$TCR!="NOT IN",c('seurat_clusters','TCR')]))
    #tmp_a<-as.data.frame(table(tmp_rds@meta.data[tmp_rds@meta.data$TCR!="NOT IN",c('seurat_clusters','TCR')]))
    tmp_a$seurat_clusters<-as.numeric(tmp_a$seurat_clusters)
    tmp_a$TCR<-i
    setwd(outdir)
    write.table(tmp_a,paste("top",which(top10_cdr3==i),"_cdr3.Freq.xls",sep=""),quote=F,sep="\t",row.names=F)
    ###样品分布
    tmp_a<-as.data.frame(table(tmp_rds@meta.data[tmp_rds@meta.data$TCR!="NOT IN",c('stim','TCR')]))
    #tmp_a$seurat_clusters<-as.numeric(tmp_a$seurat_clusters)
    tmp_a$TCR<-i
    write.table(tmp_a,paste("top",which(top10_cdr3==i),"_cdr3.sample.Freq.xls",sep=""),quote=F,sep="\t",row.names=F)
  }
}

#' Title  circlize_plot                     VDJ gene usage circlize_plot
#' @param rds                               seurat object
#' @param all_samples                       Plot vj gene usage for every samples
#' @param prefix                            Output file name prefix
#' @param outdir                            Outdir for outfile
#' @param top_n                             top_n gene for plot,default 10
#' @param rid.col                           gene color for plot,default null

#'
#' @return
#' @export
#' @importFrom   circlize  chordDiagram
#' @importFrom   circlize  circos.par
#' @importFrom   dplyr     arrange
#' @importFrom   Seurat    DimPlot
#' @importFrom   cowplot   plot_grid
#' @importFrom   ggplot2  theme
#' @examples     circlize_plot(rds,all_samples,prefix,outdir)
#'
circlize_plot<-function(rds,all_samples,prefix,outdir,top_n=10,grid.col=NULL){
  meta<-rds@meta.data[rds@meta.data$cdr3!="None",]
  vdj_tra_trb<-unique(meta[,c('cdr3','TRA','TRB')])
  cdr3_list<-c()
  for( i in 1:length(all_samples)){
    tmp_cdr3<-as.data.frame(table(meta[meta$stim==all_samples[i],'cdr3']))
    tmp_cdr3<-dplyr::arrange(as.data.frame(table(meta[meta$stim==all_samples[i],'cdr3'])), desc(Freq))
    colnames(tmp_cdr3)[1]<-"cdr3"
    tmp_cdr3<-merge(tmp_cdr3,vdj_tra_trb)
    tmp_cdr3<-dplyr::arrange(tmp_cdr3, desc(Freq))
    ####输出每个样品的TCR
    tmp_cdr4<-tmp_cdr3
    tmp_cdr4$samples<-rep(all_samples[i],nrow(tmp_cdr4))
    print(head(tmp_cdr4))
    #all_samples_tcr<<-rbind(all_samples_tcr,tmp_cdr4)
    cdr3_list<-c(cdr3_list,tmp_cdr3[1:as.numeric(top_n),c('TRA','TRB','Freq')]$TRA,tmp_cdr3[1:as.numeric(top_n),c('TRA','TRB','Freq')]$TRB)
    pdf(paste(outdir,paste(prefix,all_samples[i],"circos.pdf",sep="_"),sep='/'),w=25,h=20)#
    par(cex = 2,las=3)
    circos.par("track.height" = 0.3)
    chordDiagram(tmp_cdr3[1:as.numeric(top_n),c('TRA','TRB','Freq')],grid.col=grid.col,annotationTrack = "grid", preAllocateTracks = 1,big.gap = 10,annotationTrackHeight = c(0.1, 0.05)) #chordDiagram(mat, order = c("S2", "S1", "S3", "E4", "E1", "E5", "E2", "E6", "E3"))
    circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
      xlim = get.cell.meta.data("xlim")
      ylim = get.cell.meta.data("ylim")
      sector.name = get.cell.meta.data("sector.index")
      circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
      circos.axis(h = "top", labels.cex = 0.5, major.tick.percentage = 0.2, sector.index = sector.name, track.index = 2)
    }, bg.border = NA)
    title(all_samples[i],cex.lab=20)
    dev.off()
    write.table(tmp_cdr3,paste(outdir,paste(prefix,all_samples[i],"circos.xls",sep="_"),sep='/'),sep="\t",quote=F,row.names=F)
    tmp_cdr3<-tmp_cdr3[,c('cdr3','Freq')]
    colnames(tmp_cdr3)<-c('cdr3',all_samples[i])
  }
  return (cdr3_list)
}

#' Title  circlize_plot                     VDJ gene usage circlize_plot
#' @param immune.combined                   seurat object
#' @param all_samples                       Plot vj gene usage for every samples
#' @param prefix                            Output file name prefix
#' @param outdir                            Outdir for outfile
#' @param rid.col                           gene color for plot,default null

#'
#' @return
#' @export
#' @importFrom   circlize       chordDiagram
#' @importFrom   circlize       circos.par
#' @importFrom   dplyr          arrange
#' @importFrom   circlize       rand_color
#' @importFrom   cowplot        plot_grid
#' @importFrom   ggplot2  theme
#' @examples     circ_plot(immune.combined,outdir)
#'
circ_plot<-function(immune.combined,outdir,prefix='gene_using',top_n=10,all_samples=NULL,grid.col=NULL){
  library(circlize)
  mkdirs(outdir,'circ_plot')
  setwd(paste(outdir,'circ_plot',sep='/'))
  if (is.null(all_samples)){all_samples=unique(immune.combined@meta.data$stim)}
  cdr3_list<-circlize_plot(immune.combined,all_samples,prefix,paste(outdir,'circ_plot',sep='/'),top_n=10)
  grid.col<-rand_color(length(cdr3_list))
  names(grid.col)<-cdr3_list
  cdr3_list<-circlize_plot(immune.combined,all_samples,prefix,paste(outdir,'circ_plot',sep='/'),top_n=10,grid.col=grid.col)

}

#' Title  circlize_plot                     VDJ gene usage circlize_plot
#' @param tcr                               tcr freq data.frame
#' @param outfile                           Output heatmap file name
#' @param w_h                               Output file w and h
#' @param topn                              Top n cdr3 for heatmap plot
#' @param title                             title for plot file

#'
#' @return
#' @export
#' @importFrom   pheatmap       pheatmap
#' @importFrom   ggplot2  theme
#' @examples     tcr_heatmapplot(tcr)
#'
tcr_heatmapplot<-function(tcr,outfile='CDR3_counts_samples.TCRnames.pdf',w_h=c(12,10),topn=30,title='CDR3 counts'){
  pdf(outfile,w=w_h[1],h=w_h[2])
  library(pheatmap)
  tmp_tcr<-tcr[,2:(ncol(tcr)-1)]
  rownames(tmp_tcr)<-tcr[,1]
  p<-pheatmap(log(tmp_tcr[1:topn,]+1,2),show_rownames=T,show_colnames=T,cluster_cols = T,cluster_rows= T,treeheight_row=2,border=FALSE,main=title,number_color = "black",cellwidth=40,fontsize_row=15,cutree_row = 3,fontsize = 18,angle_col=45)
  print(p)
  dev.off()
}


